C  seis88 - Rays tracing tool
C  Copyright (C) 2009 Ivan Psencik <ip@ig.cas.cz>
C 
C  This program is free software: you can redistribute it and/or modify
C  it under the terms of the GNU General Public License as published by
C  the Free Software Foundation, either version 3 of the License, or
C  (at your option) any later version.
C 
C  This program is distributed in the hope that it will be useful,
C  but WITHOUT ANY WARRANTY; without even the implied warranty of
C  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C  GNU General Public License for more details.
C 
C  You should have received a copy of the GNU General Public License
C  along with this program.  If not, see <http://www.gnu.org/licenses/>.


C
C     P R O G R A M   P O L A R P L O T
C     *********************************
C
C     PROGRAM POLARPLOT IS DESIGNED FOR THE PLOTTING OF POLAR
C     DIAGRAMS FROM THE LOGICAL UNITS LU3 GENERATED BY SYNTPL
C
C     ************************************************************
C
      DIMENSION SEISA(3001),SEISB(3001),IEP(100),TEXT(17),TXT(17)
C
C
      IRUN=0
      mode=0
      call serv(mode,lin,lou,llu3a,llu3b,ldum)
      if(mode.eq.0)lin=5
      if(mode.eq.0)lou=6
      if(mode.eq.0)CALL PLOTS(ldum1,ldum2,7)
      CALL PLOT(5.,5.,-3)
    1 READ(lin,101)LU3A,LU3B,ISHIFT,IPRINT
      WRITE(lou,101)LU3A,LU3B,ISHIFT,IPRINT
      if(mode.eq.1)lu3a=llu3a
      if(mode.eq.1)lu3b=llu3b
C
C
      IF(ISHIFT.EQ.0)ISHIFT=10
      SHIFT=FLOAT(ISHIFT)
    2 READ(lin,101)MCONT,MEPIC,NTICX,NTEXT
      WRITE(lou,101)MCONT,MEPIC,NTICX,NTEXT
      IF(MCONT.EQ.0)GO TO 99
      IF(MCONT.EQ.(-1))REWIND LU3a
      IF(MCONT.EQ.(-1))REWIND LU3b
      IF(MCONT.EQ.(-1))GO TO 1
      IF(MEPIC.EQ.0)GO TO 3
      READ(lin,101)NEPIC,(IEP(I),I=1,NEPIC)
      WRITE(lou,101)NEPIC,(IEP(I),I=1,NEPIC)
    3 CONTINUE
      READ(lin,102)XLEN,DTICX,SC,TSTART,TFIN,AMP,B1
      WRITE(lou,102)XLEN,DTICX,SC,TSTART,TFIN,AMP,B1
      IF(ABS(SC).LT..00001)SC=1.
      IF(ABS(B1).LT..00001)B1=1.
      READ(lin,100)TXT
      WRITE(lou,100)TXT
      REWIND LU3A
      READ(LU3A,100)TEXT
      WRITE(lou,100)TEXT
      READ(LU3A,105)MDISTA,MRED,MCOMPA,itpr,VRED,RSTEP,XSOUR,DT
      WRITE(lou,105)MDISTA,MRED,MCOMPA,itpr,VRED,RSTEP,XSOUR,DT
      READ(LU3A,104)XMXA,AMAXIM
      WRITE(lou,104)XMXA,AMAXIM
      REWIND LU3B
      READ(LU3B,100)TEXT
      WRITE(lou,100)TEXT
      READ(LU3B,105)MDISTB,MRED,MCOMPB,itpr,VRED,RSTEP,XSOUR,DT
      WRITE(lou,105)MDISTB,MRED,MCOMPB,itpr,VRED,RSTEP,XSOUR,DT
      READ(LU3B,104)XMXB,BMAXIM
      WRITE(lou,104)XMXB,BMAXIM
      IF(MDISTA.NE.MDISTB)WRITE(lou,108)
      IF(MDISTA.NE.MDISTB)GO TO 99
      SMAXIM=AMAXIM
      IF(SMAXIM.LT.BMAXIM)SMAXIM=BMAXIM
C
C     LOOP FOR THE RECEIVER POSITIONS
C
      DO 10 I=1,MDISTA
      READ(LU3A,110)XX,SMAXIA,TMINA,NPTSA
      READ(LU3A,109)(SEISA(M),M=1,NPTSA)
      READ(LU3B,110)XX,SMAXIB,TMINB,NPTSB
      READ(LU3B,109)(SEISB(M),M=1,NPTSB)
      IF(MEPIC.EQ.0)GO TO 5
      DO 6 J=1,NEPIC
      IF(I.EQ.IEP(J))GO TO 5
    6 CONTINUE
      GO TO 10
    5 SAUX=SMAXIA/999.
      DO 22 M=1,NPTSA
   22 SEISA(M)=SEISA(M)*SAUX
      SAUX=SMAXIB/999.
      DO 23 M=1,NPTSB
   23 SEISB(M)=SEISB(M)*SAUX
C
C     PLOT OF FRAME
      XMER=.5*XLEN
      DDX=XMER
      IF(IRUN.NE.0)CALL PLOT(XLEN+SHIFT,0.,-3)
      IRUN=1
C*..................
      CALL COLOR(14)
C*..................
      IF(NTEXT.NE.0)CALL BORDER(XLEN,DTICX,XLEN,DTICX,SC,TXT,0,-1.,
     11.,-1.,1.,NTICX,NTICX,1,1)
      IF(NTEXT.EQ.0)CALL BORDER(XLEN,DTICX,XLEN,DTICX,SC,TEXT,0,
     1-1.,1.,-1.,1.,NTICX,NTICX,1,1)
      ELM=.45*SC
      T=.5*(XLEN-6.*ELM)
      IF(MCOMPA.EQ.1)
     1CALL SYMBOL(T,-1.6*SC,ELM,'RADIAL',0.,6)
      T=T-ELM
      IF(MCOMPA.EQ.0)
     1CALL SYMBOL(T,-1.6*SC,ELM,'VERTICAL',0.,8)
      T=T-ELM
      IF(MCOMPA.EQ.2)
     1CALL SYMBOL(T,-1.6*SC,ELM,'TRANSVERSE',0.,10)
      U=-(1.6+.4)*SC
      IF(MCOMPB.EQ.2)
     1CALL SYMBOL(U,T,ELM,'TRANSVERSE',90.,10)
      T=T+ELM
      IF(MCOMPB.EQ.0)
     1CALL SYMBOL(U,T,ELM,'VERTICAL',90.,8)
      T=T+ELM
      IF(MCOMPB.EQ.1)
     1CALL SYMBOL(U,T,ELM,'RADIAL',90.,6)
      CALL PLOT(0.,0.,3)
      if(itpr.ne.22)CALL SYMBOL(ELM,XLEN+SC,ELM,'3HX= ',0.,3)
      if(itpr.eq.22)CALL SYMBOL(ELM,XLEN+SC,ELM,'3HZ= ',0.,3)
      CALL NUMBER(4*ELM,XLEN+SC,ELM,XX,0.,2)
      CALL SYMBOL(10.*ELM,XLEN+SC,ELM,'KM, ',0.,4)
      CALL SYMBOL(14.*ELM,XLEN+SC,.5*ELM,'T1,T2= ',0.,7)
      CALL NUMBER(17.5*ELM,XLEN+SC,.5*ELM,TSTART,0.,2)
      CALL NUMBER(20.*ELM,XLEN+SC,.5*ELM,TFIN,0.,2)
C
C
      SMAXI=SMAXIA
      IF(SMAXI.LT.SMAXIB)SMAXI=SMAXIB
      TMIN=TMINA
      IF(TMIN.GT.TMINB)TMIN=TMINB
      TMAXA=TMINA+(NPTSA-1)*DT
      TMAXB=TMINB+(NPTSB-1)*DT
      TMAX=TMAXA
      IF(TMAX.LT.TMAXB)TMAX=TMAXB
      if(iprint.ge.2)WRITE(lou,102)TMIN,TMAX,SMAXI,SMAXIM
C
      IF(SMAXI.LT.0.000001)GO TO 7
      IF(ABS(AMP).LT.0.00001)FACTOR=B1*DDX/SMAXI
      IF(ABS(AMP).LT.0.00001)GO TO 8
      IF(AMP.LT.(-0.00001))FACTOR=B1*DDX/SMAXIM
      IF(AMP.GT.0.00001)FACTOR=B1
      GO TO 8
    7 FACTOR=0.
    8 CONTINUE
      SFMAX=FACTOR*SMAXI
      SF1=.003*SFMAX
      IF(IPRINT.ge.1)WRITE(lou,103)XX,SMAXI,FACTOR,SFMAX
C*.................
      CALL COLOR(4)
C*.................
      K=0
      IA=0
      IB=0
      XNEW=0.
      YNEW=0.
      HLEN=.5*XLEN
      TST=TSTART
      TEND=TFIN
      IF(TST.LT.TMIN)TST=TMIN
      IF(TEND.GT.TMAX)TEND=TMAX
      IF(TST.LT.TMINA)XNEW=0.
      IF(TST.LT.TMINA)GO TO 14
      IA=(TST-TMINA)/DT+1
      T=TMINA+DT*FLOAT(IA-1)
      TM=T
      AMPL=SEISA(IA)+(SEISA(IA+1)-SEISA(IA))*(TST-T)/DT
      XNEW=FACTOR*AMPL
      IF(ABS(XNEW).GT.HLEN)GO TO 15
   14 IF(TST.LT.TMINB)YNEW=0.
      IF(TST.LT.TMINB)GO TO 12
      IB=(TST-TMINB)/DT+1
      T=TMINB+DT*FLOAT(IB-1)
      TM=T
      BMPL=SEISB(IB)+(SEISB(IB+1)-SEISB(IB))*(TST-T)/DT
      YNEW=FACTOR*BMPL
      IF(ABS(YNEW).GT.HLEN)GO TO 15
   12 CONTINUE
      XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      if(iprint.ge.2)WRITE(lou,102)XNEW,YNEW,AMPL,BMPL
      CALL PLOT(XNEW,YNEW,3)
   15 CONTINUE
      IF(ABS(T-TMINA).LT..001)IA=1
      IF(ABS(T-TMINB).LT..001)IB=1
      IF(IA.NE.0)IA=IA+1
      IF(IB.NE.0)IB=IB+1
      XNEW=0.
      IF(IA.GT.0.AND.IA.LE.NPTSA)XNEW=FACTOR*SEISA(IA)
      YNEW=0.
      IF(IB.GT.0.AND.IB.LE.NPTSB)YNEW=FACTOR*SEISB(IB)
      IF(ABS(XNEW).GT.HLEN)GO TO 15
      IF(ABS(YNEW).GT.HLEN)GO TO 15
      XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      K=K+1
      T=TM+K*DT
      IF(T.GT.TEND)GO TO 13
      if(iprint.ge.2)WRITE(lou,106)K,IA,IB,T,XNEW,YNEW
      CALL PLOT(XNEW,YNEW,2)
      GO TO 15
   13 XNEW=0.
      IF(T.GT.TMAXA)GO TO 11
      AMPL=SEISA(IA-1)+(SEISA(IA)-SEISA(IA-1))*(TMAXA-T)/DT
      XNEW=FACTOR*AMPL
      IF(ABS(XNEW).GT.HLEN)GO TO 10
   11 YNEW=0.
      IF(T.GT.TMAXB)GO TO 9
      BMPL=SEISB(IB-1)+(SEISB(IB)-SEISB(IB-1))*(TMAXB-T)/DT
      YNEW=FACTOR*BMPL
      IF(ABS(YNEW).GT.HLEN)GO TO 10
    9 XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      CALL PLOT(XNEW,YNEW,2)
   10 CONTINUE
C
C     END OF THE LOOP FOR RECEIVER POSITIONS
C
      GO TO 2
C
C
  100 FORMAT(17A4)
  101 FORMAT(16I5)
  102 FORMAT(8F10.5)
  103 FORMAT(2X,4E15.5)
  104 FORMAT(22X,F10.5,9X,E15.9)
  105 FORMAT(4I5,4F10.5)
  106 FORMAT(3I5,4F10.5)
  108 FORMAT(/1X,'DIFFERENT SELECTION OF RANGES ON THE AXES,
     1 COMPUTATION TERMINATED'//)
  109 FORMAT(20F4.0)
  110 FORMAT(F10.5,E15.8,F10.5,I5)
   99 CALL PLOT(0.,0.,999)
C
C
      STOP
      END
